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The study of incompressible magnetohydrodynamic (MHD) turbulence gives useful insights on many astrophys- 
ical problems. We describe a pseudo-spectral MHD code suitable for the study of incompressible turbulence. We 
review our recent works on direct three-dimensional numerical simulations for MHD turbulence in a periodic box. 
In those works, we use a pseudo-spectral code to solve the incompressible MHD equations. We first discuss the 
structure and properties of turbulence as functions of scale. The results are consistent with the scaling law recently 
proposed by Goldreich & Sridhar. The scaling law is based on the concept of scale-dependent isotropy: smaller ed- 
dies are more elongated than larger ones along magnetic field lines. This scaling law substantially changes our views 
on MHD turbulence. For example, as noted by Lazarian & Vishniac, the scaling law can provide a fast reconnec- 
tion rate. We further discuss how the study of incompressible MHD turbulence can help us to understand physical 
processes in interstellar medium (ISM) by considering imbalanced cascade and viscous damped turbulence. 
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I. INTRODUCTION 

Astrophysical flows are complicated and dynamic. 
For example, molecular clouds are clumpy over a vast 
range of scales and observations of velocity spectral line 
broadening indicate that dynamic pressure {p-v^/2) 
usually dominates the thermal pressure (nkT). Fur- 
thermore, velocity spectral lines show multiple compo- 
nents. These facts suggest that cloud internal velocity 
possesses disordered components. In molecular clouds, 
the characteristic Reynolds number is estimated to be 
larger than ~ 10^°, which is consistent with the turbu- 
lent state in the clouds. In general, when the Reynolds 
number (= LV/v; L=characteristic size of the system, 
V=velocity dispersion, and i^=viscosity) is larger than 
an order of 10 - 100, the system becomes unstable and 
turbulent. 

In laboratory systems, fluids are usually incompress- 
ible and unmagnetized. Kolmogorov phenomenology 
provides an excellent insight for such systems. Sup- 
pose that we 'disturb' the fluid at a scale L. We call 
this scale as energy injection scale or largest energy con- 
taining eddy scale. Then the energy injected to the 
scale L cascades to progressively smaller and smaller 
scales. Ultimately, the energy will reach the molecular 
dissipation scale Id and it will be lost there. The scales 
between L and Id is called the inertial range. Suppose a 
scale I lies in the inertial range. Let the characteristic 
velocity associated with the scale be vi. Kolmogorov 
theory states that the kinetic energy (vf) is transferred 
to one-level smaller scale within one eddy turnover time 
(l/vi). It is natural that the cascade rate {vf /{l/vi)) be 
scale-independent. Therefore, we have 



Vl cx 



l^/\ 



(1) 



Note that E(k)dk is the amount of energy between 
the wavenumber k and k + dk. When E{k) follows 
a power law, kE(k) is the energy near the wavenumber 
k ~ l/l. Since vf represents a similar energy, we have 
vf K. kE{k). Therefore, equation (P becomes 



E{k) cx k-^/^. 



(2) 



(One-dimensional) Energy spectrum E[k) is one of 
the most important quantities in turbulence theories. 



This is the well-known Kolmogorov spectrum. Astro- 
physical fluids are compressible and magnetized. Nev- 
ertheless, many astrophysical quantities follow this Kol- 
mogorov spectrum. It is important to note that Kol- 
mogorov theory predicts: 1. the velocity dispersion 
(^ Vl) decreases with scales as vi oc V-^'^; 2. the slope 
of the energy spectrum is —5/3. 

Earlier works on ISM focused on the velocity disper- 
sions and derived a qualitative agreement with the Kol- 
mogorov theory. Larson (1981) and Myers (1983) found 
that velocity dispersion decreases as the size of the 
clumps decreases, which may support turbulent model 
of ISM. Scalo (1984) showed that two-point velocity 
correlation increases as the separation decreases, which 
is also qualitatively consistent with turbulent model. 
All the above mentioned studies are for scales larger 
than or similar to ~ 1 pc. 

On the other hand, other studies address the issue 
of power spectra. Interstellar scintillation observations 
indicate electron density spectrum follows a power law 
over 7 decades of length scales (see Armstrong et al. 
1995). The slope of the spectrum is very close to —5/3 
for lO^m - lO^^m. The spectrum is somewhat uncertain 
for scales larger than 10^'^m (~ O.Olpc). However, re- 
cently, Lazarian & Pogosyan (2000) and Starnimirovic 
& Lazarian (2001) showed the existence of the —5/3 
velocity power spectrum over pc-scales in HI. Solar- 
wind observations provide in-situ measurements of the 
power spectrum of magnetic fluctuations. Leamon et 
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al. (1998) obtained a slope of —1.7. 

Therefore it is evident that turbulence is prevalent in 
astrophysical fluids. Understanding the nature of such 
turbulence is important because it affects many astro- 
physical processes. For example, star formation rate 
depends on the turbulence decay time-scale. Heat and 
cosmic ray transport also depend on the statistics of 
turbulence. Recently anisotropic MHD Turbulence has 
also been evoked to provide fast magnetic reconnection 
rate (Lazarian & Vishniac 1999). 

One approach to studying the ISM is to perform 
time-dependent numerical simulations to model the 
ISM, including as many of the interacting phenomena 
as practical. Of course, the physics included in such 
models must necessarily be highly simplified, and it is 
difficult to determine which features of the final model 
result from which physical assumptions (or initial con- 
ditions). Our approach is to use simplified numerical 
simulations to study the influences of various physical 
phenomena in isolation. We want to obtain a physical 
feeling for the general cfi^ects that each phenomenon has 
on the nature of the ISM. Here we review our previous 
works on incompressible MHD turbulence. It is obvi- 
ous that the real ISM is compressible, but we want to 
separate the effects of magnetic turbulence from those 
involving compression. In later works we will include 
compression for comparison with the present models, 
thereby isolating its importance directly. 

As we said earlier, we will numerically study incom- 
pressible MHD turbulence. We believe Fourier spectral 
method is one of the best numerical methods for such 
calculations. In II, we show how the spectral method 
works. In III, we introduce our recent works done with 
the spectral code: anisotropy, decay speed of MHD tur- 
bulence, and magnetic structures below the viscous cut- 
off. In IV, we give conclusion. 

II. NUMERICAL METHOD 

(a) PSEUDO-SPECTRAL METHOD 

For simplicity, let us flrst consider incompressible 
unmagnetized fluid. Then, the governing equation is 
the Navier-Stocks equation: 



dt 



-V • (vv) + i/V^v + VP, 



(3) 



where v is velocity, v is viscosity, and P is gas pressure. 
Since we are considering an incompressible fluid, we can 
assume that p= 1- 

The Fourier spectral method is highly suitable for in- 
compressible flows. Simply put, the idea of the method 
is simple: we transform all physical space variables into 
Fourier (or wave- vector) space and solve the differential 
equation there. Hereafter, we assume that we are solv- 
ing the equation on a uniform rectangular grid. Some 
advantage of the method is: 

1. It shows fast convergence for smooth functions. 

2. It is easy to implement divergence-free conditions. 



3. It is easy to implement periodic boundary condition. 

4. Fourier space analyses (e.g. calculation of energy 
spectrum) are natural consequences. 

We note that the gradient operator becomes ik. and 
that multiplicatio n be comes convolution sum in Fourier 
space. Here i = -\/^. Using these facts, we obtain the 
Navier-Stocks equation in Fourier space: 



9v(k) 

dt 



-ik • (w)(k) - iukH{k) + ikP(k), (4) 



where 



(vv)(k)= / dkv(p)v(q). (5) 

J p+q=k 

The incompressibility condition becomes 



ik ■ v(k) = 0, 



(6) 



which means that the velocity vector stays perpendic- 
ular to the wave vector k. Since the pressure term in 
equation (4) is parallel to k, we can drop the term. On 
the other hand, the nonlinear term (i.e the term con- 
taining vv) may have both parallel and perpendicular 
components. Therefore, we need to remove the compo- 
nent parallel to k by applying an appropriate projection 
operator. 

The calculation of the nonlinear term is the most 
time consuming part of the Fourier spectral methods. 
For simplicity, let us consider the following one dimen- 
sional convolution sum: 



dk HpMq)- 



p+q=k 



(7) 



Then the discrete representation of the convolution sum 
is 

p+g=k 

The pseudo-spectral method approximates the convo- 
lution sum as 



p+q=k 



(9) 



where F and denote forward/inverse Fourier trans- 
form. The cost of such calculation is 0{N log N), where 
iV is the number of grid points. 

Once we obtain the convolution sum, we can use any 

techniques for ordinary differential equations to solve 
the time evolution of the Navier-Stocks equation. 

(b) ACTUAL SETUP 

In actual calculations, we have adopted a pseu- 
dospectral code and calculated the time evolution of 

incompressible magnetic turbulence subject to a ran- 
dom driving force per unit mass: 

^ = (Vxv)xv-(VxB)xB-Fi/V^v-|-f-FVP', (10) 
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V X (v X B) + r?V^B, (11) 

V-v = V-B = 0, (12) 

where f is a random driving force, P' = P/p + v ■ v/2, 

V is the velocity, and B is magnetic field divided by 
{AnpY^^. In this representation, v can be viewed as 
the velocity measured in units of the r.m.s. velocity, v, 
of the system and B as the Alfven speed in the same 
units. The time t is in units of the large eddy turnover 
time (~ L/v) and the length in units of L, the inverse 
wavenumber of the fundamental box mode. In this sys- 
tem of units, the viscosity and magnetic diffusivity 
T] are the inverse of the kinetic and magnetic Reynolds 
numbers respectively. The turbulence is driven by 21 
random forcing components in Fourier space. The peak 
of energy injection occurs at fc ~ 2.5. The amplitudes 
of the forcing components are tuned to ensure u « 1 
We use exactly the same forcing terms for all simula- 
tions. We use up to 256^ collocation points. We use an 
integration factor technique for kinetic and magnetic 
dissipation terms and a leap-frog method for nonlinear 
terms. At t = 0, the magnetic field has only its uni- 
form component and the velocity field is restricted to 
the range 2 < fc < 4 in wavevector space. 

For some runs, we use hyperviscosity and hypcrdif- 
fusivity for the dissipation terms. The power of hyper- 
viscosity is set to either 2 or 8, so that the dissipation 
term in the above equation is replaced with 



100 



- t^/.(v2)V, 



(13) 



where /i = 2 or 8. A similar expression is used for the 
magnetic dissipation term. In subsections III (a) and 
Ill(b), we use unit magnetic Prandtl number {v = rj) 

and the notation NY-_BoZ. In subsection III(c), we use 
physical viscosity and physical- or hypcr-difFusion and 
the notation NXY-SoZ. Here N = 256, 144 refers to the 
number of grid points in each spatial direction; X = P 
refers to physical viscosity; Y = P, H2, H8 refers to the 
form of diffusion (physical-, the 2nd power hyper-, or 
the 8th power hyper-diffusion, respectively); Z=0.5 or 
1 refers to the strength of the external magnetic field 
in Alfven velocity unit. 

Diagnostics for our code can be found in Cho and 
Vishniac (2000b). 

III. PROPERTIES OF MHD TURBULENCE 

In this section, we discuss the properties of MHD 
turbulence in the presence of a strong uniform back- 
groimd field. 

(a) ANISOTROPY 

Historically hydrodynamic turbulence in an incom- 
pressible fluid was successfully described by the eddy 
cascade (Kolmogorov 1941), but MHD turbulence was 
first modeled by wave turbulence (Iroshnikov 1963, 
Kraichnan 1965). This theory assumes isotropy of the 
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v: 256H-Bj,1 
b: 256H-Bj,1 
v; 256H-B„0.5 
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b: 256P-B 1 
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Semi-minor Axis (~1/I<j^) 



Fig. 1. — The relation between fc|| and k_\_. (Cho & 
Vishniac 2000a). 



energy cascade in Fourier space. However, the mag- 
netic field defines a local symmetry axis since it is easy 
to mix field lines in directions perpendicular to the local 
B and much more difficult to bend them. In a turbulent 
medium, the kinetic energy associated with large scale 
motions is greater than that of small scales. However, 
the strength of the local mean magnetic field is almost 
the same on all scales. Therefore, it becomes relatively 
difficult to bend magnetic field lines as we consider 
smaller scales, leading to more pronounced anisotropy. 
A self-consistent model of MHD turbulence which in- 
corporates this concept of scale dependent anisotropy 
was introduced by Goldreich & Sridhar (1995; here- 
inafter GS95). 

The major predictions of the GS95 model are as fol- 
lows. First, anisotropy is scale dependent. Roughly 
speaking, the semi-major axis (l/fc||) and the semi- 
minor axis {\/k±^ of an eddy satisfy the relation k\\ oc 

A;^^"^, which means smaller eddies are relatively more 
elongated than larger ones. Here k± and fcy are 
wave numbers perpendicular and parallel to the back- 
ground field. Second, the model predicts that the one- 
dimensional energy spectrum is Kolmogorov-typc if ex- 
pressed in terms of perpendicular wavenumbers, i.e. 
E{ki_) oc k~^_^'^. 

Numerical simulations by Cho & Vishniac (2000a) 
and Maron & Goldreich (2001) have mostly supported 
the GS95 model and helped to extend it. Both analyses 
stressed the point that scale dependent anisotropy can 
be measured only in local coordinate frames which are 
aligned with the locally av(a'agc;d magnetic field direc- 
tion. Cho & Vishniac (2000a) calculated the structure 
functions of velocity and magnetic field in the local 
frames, and found that the contours of the structure 
functions do show scale dependent anisotropy, consis- 
tent with the predictions of the GS95 model. Fig 1. 
shows that the semi-major axis (l/fcy) is proportional 



4 



CHO 



to the 2/3 power of the semi-minor axis {l/k±), imply- 



ing that kn 



u2/3 



One dimensional energy spectrum 



follows Kolmogorov spectrum: E{k) oc k (see Cho 
& Vishniac 2000a). 

(b) DECAY OF MHD TURBULENCE 

Turbulence plays a critical role in molecular cloud 

support and star formation and the issue of the time 
scale of turbulent decay is vital for understanding these 
processes. 

If MHD turbulence decays quickly then serious prob- 
lems face researchers attempting to explain important 
observational facts, i.e. turbulent motions scon within 
molecular clouds without star formation (sec Myers 
1999) and rates of star formation (Mckee 2000). Ear- 
lier studies attributed the rapid decay of turbulence 
to compressibility effects (Mac Low 1999). Cho et al. 
(2001), as well as earlier ones (Cho & Vishniac 2000a, 
Maron & Goldreich 2001), shows that turbulence de- 
cays rapidly even in the incompressible limit. This can 
be understood in the framework of GS95 model where 
mixing motions perpendicular to magnetic field lines 
form hydrodynamic-typc eddies. Such eddies, as in 
hydrodynamic turbulence, decay in one eddy turnover 
time. 

Below we consider the effect of imbalance on the 
turbulence decay time scale. MHD turbulence can be 
described by opposite-traveling wave packets. 'Imbal- 
ance' means that wave packets traveling in one direc- 
tion have significantly larger amplitudes than the other. 
In astronomy, many energy sources are localized. For 
example, SN explosions and OB winds are typical point 
energy sources. Furthermore, astrophysical jets from 
YSOs are believed to be highly collimated. With these 
localized energy sources, it is natural to think that in- 
terstellar turbulence is typically imbalanccd. 

In this subsection, we explicitly relate the degree of 
imbalance and the decay time scale of turbulence in the 
presence of a strong uniform background field. 

In Figure 2 we demonstrate that an imbalanccd cas- 
cade does extend the lifetime of MHD turbulence. We 
use the run 144H-i?ol to investigate the decay time 
scale. We ran the simulation up to t=75 with non-zero 
driving forces. Then at t=75, we turned off the driving 
forces and let the turbulence decay. At t=75, there is a 
slight imbalance between upward and downward mov- 
ing components {E+ = 0.499 and E_ = 0.40). This re- 
sults from a natural fluctuation in the simulation. The 
case of (i^-i-)to = 80%(£J_)to corresponds to the simula- 
tion that starts off from this initial imbalance. In other 
cases, wc either increase or decrease the energy of "l 
components and, by turning off the forcing terms, let 
the turbulence decay. We can clearly observe that im- 
balanccd turbulence extends the decay time scale sub- 
stantially. Note that we normalized the initial energy 
to 1. The; y-axis is the total (=up + down) energy. 

In this section, we found that turbulence decay time 




Fig. 2. — Decay of imbalanccd turbulence. Imbal- 
anccd cascade can extend decay time. Run 144H-Bol. 
(Cho et al. 2001) 



can be slow. If wc consider the interstellar medium at 
larger scales, it is definitely compressible, has a whole 
range of energy injection/dissipation scales (see Scalo 
1987), and the relative role of vortical versus com- 
pressible motions being unclear. Nevertheless, we be- 
lieve that our simplified treatment may still elucidate 
some of the basic processes. To what extend this claim 
can be justified will be clear when we compare com- 
pressible and incompressible results. However, if we 
accc^pt that fast and slow MHD modes are subjected 
to fast coUisionless damping (see Minter & Spangler 
1997) the remaining modes are incompressible Alfven 
modes. Those should be well described by our model 
when turbulence is supersonic but sub-Alfvenic. Our 
preliminary results (Cho & Lazarian 2002) show that 
the coupling of the modes is marginal even in compress- 
ible regime. 

(c) STRUCTURES BELOW VISCOUS CUT- 
OFF 

In hydrodynamic turbulence viscosity sets a mini- 
mal scale for motion, with an exponential suppression 
of motion on smaller scales. Below the viscous cutoff 
the kinetic energy contained in a wavcnumber band is 
dissipated at that scale, instead of being transferred 
to smaller scales. This means the end of the hydro- 
dynamic cascade, but in MHD turbulence this is not 
the end of magnetic structure evolution. For viscosity 
much larger than resistivity, v ^ j], there will be a 
broad range of scales where viscosity is important but 
resistivity is not. On these scales magnetic field struc- 
tures will be created through a combination of large 
scale shear and the small scale motions generated by 
magnetic tension. As a result, we expect a power-law 
tail in the energy distribution, rather than an exponen- 
tial cutoff. To our best knowledge, this is a completely 
new regime of MHD turbulence. Ambipolar diffusion 
damps turbulent motions in ISM. Therefore, the study 
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that magnetic fields have stochastic structures below 
the ambipolar diffusion scale in ISM and, therefore, af- 
fect many astrophysical processes, such as cosmic ray 
transports, magnetic reconnection, and thermal diffu- 
sion. 
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of such a viscous damped MHD turbulence is particu- 
larly important for many astrophysical processes below 
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In Cho, Lazarian, & Vishniac (2002), we numeri- 
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netic energy spectrum below the viscous damping scale. 
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different values of magnetic diffusion: physical diffu- 
sion, hyper-diffusion with order 2, and hyper-diffusion 
with order 8. We will present a theoretical model for 
this regime in an upcoming paper (Lazarian, Vishniac, 
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corresponds to the energy injection scale. Second, for 
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in general sense. Third, magnetic and kinetic spectra 
decouple at /c ^ 6. Forth, for 20 < /c, it appears that 
a new damped-scale inertial range emerges. The emer- 
gence of the new inertial range is conspicuous in the 
order 8 hyperdiffusion simulation. 

IV. SUMMARY 

We described a pseudo-spectral MHD code and demon- 
strated its usefulness in astrophysics through our recent 
works. We showed that the rate at which MHD turbu- 
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ance between wave packets traveling in opposite direc- 
tions. A substantial degree of imbalance can substan- 
tially extend the decay time scale of the MHD turbu- 
lence, which might be useful for the study of star forma- 
tion. We also showed that magnetic fields can have rich 
structures below the viscous cut-off scale, which implies 
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